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Abstract 


Crystallization of the lunar magma ocean creates a chemically stratified Moon consisting of an anorthositic crust and 

O’vstalhze^o" T , ° Verly,n ® ' he pnm “ ive lunar interior - Within magma ocean cumulates the last liquids to 
under vinp l 7 > ""*»««-"<* cumulates that contain high concentrations of incompatible radioactive elements The 

und lying olivine-orthopyrorcene cumulates are also stratified with later crystallized, denser, more Fe-rich comTsftions* 

mss hf |' S P rT r CXP eS C Chem,Cal and ,hermal consc< I uences of an internal evolution model accounting for the 
possible role of these sources of chemical buoyancy. Rayleigh-Taylor instability causes the dense ilmenite-rich cumulate 
layer and underlying Fe-rich cumulates to sink toward the center of the Moon, forming a dense lunar rnre After 
model' radl0ac,,ve heat ' n 8 within the ilmenite-rich cumulate core heats the overlying mantle, causing it to melt In this 
modd. the source region for h.gh-Ti0 2 mare basalts is a convectively mixed layer above the core-manL ^uSv whirt 
would con am small and variable amounts of admixed ilmenite and KREEP. This deep high-pressure melting as required fw 

wmfnroTSir' 6 . af,er 3 reaSOnable ,ime in,erval to «pl™ 'he onset of mare basalt volcanL if the 

Trance T a " d ‘ he chemical densi 'y gradients above the core are sufficiently high bu, wilhm 

I'fZ 7 h3t ?! gh ' h3Ve 6660 pfeSem in ,he Moon - R ^' d '<^ of details implied by particular modfoarame^ 

La ° Ver !| Urn d T" by thC h ' fih denS ' ,y ° f magma 0cean Fc - rich cumu 'ates should Concentrate high-TiO mare 

7 rCeS ;. and Pr ° f bab ^ 3 SigniflCan ' fraC,i °" ° f radioactive heali "g. 'oward the center of the Z Sis wifi hav" 
important implications for both the thermal evolution of the Moon and for mare basalt genesis. 


1. Introduction 

This study explores a model for the chemical 
evolution of the lunar interior that attempts to ex- 
plain the origin and evolution of lunar magmatism 
and possibly the existence of a lunar core [1,2], 
Inferred source depths greater than 400-500 km (cf. 
[3], [4], [5]) as well as the ages (3.9-3. 1 Ga) of high- 
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and low-Ti mare basalt volcanism [6] must reflect 
processes of lunar internal evolution. However, these 
generally accepted characteristics of lunar magmatic 
evolution do not appear to be explained by previous 
models for the evolution of the lunar interior. Many 
previous studies have focussed on geophysical (e.g., 
[7,8]) and tectonic constraints [9-13] on lunar evolu- 
tion. Among those that specifically address the solid- 
ification of the magma ocean (e.g., [14-16]), the 
study of Solomon and Longhi [14] couples cooling 
and solidification of the magma ocean with the 
larger scale lunar interior thermal evolution. Accord- 


00 1 2-821 X/95/509.50 © 1995 Elsevier Science B.V. All rights reserved 
SSDI 00 12-821 X(95)00138-7 



ing to their results, magma ocean mafic cumulates, 
the source of mare basalts required by geochemical 
constraints, never remelt. Melting of the primitive 
lunar mantle beneath the magma ocean cumulates, as 
predicted by their model, would produce magmas but 
these would not be mare basalts. The obvious con- 
clusion is that the production of mare basalts must 
involve factors not included in the Solomon and 
Longhi model. 

We accept, as a working hypothesis, that a large 
magma ocean, several hundreds of kilometers in 
thickness, segregated into a plagioclase-rich crust 
and a complimentary stratified cumulate mantle en- 
riched in Fe-Mg minerals (see [17,18]). The initial 
stratification of the magma ocean cumulate can only 
be estimated given the possible complexities associ- 
ated with the crystallization of a magma body of 
such extraordinary size (e.g., [19,20]). The solidifica- 
tion of this magma ocean could be approximated by 
models of near-equilibrium crystallization to near- 
perfect fractional crystallization, or some combina- 
tion of both processes occurring at different stages. 
In addition, cooling caused by the rapid transport of 
quenched melt to the bottom of the magma ocean 
would induce olivine crystallization, producing a 
cumulate layer approximating the composition of an 
olivine-rich bulk Moon [21]. 

There is dear evidence of fractional crystalliza- 
tion, however, at least in the final stages of solidifi- 
cation of the magma ocean. The formation of the 
anorthosite-rich crust by 4.5 Ga and the establish- 
ment of an ilmenite-bearing cumulate mantle layer 
by 4.4 Ga [6] require very effective crystal-liquid 
segregation. Calculations of the rare earth composi- 
tion of magmas parental to ferroan anorthosites sug- 
gest that the plagioclase began to segregate from the 
magma ocean only after 80-90% solidification [22]. 
Thus, fractional crystallization must have been oper- 
ative at this late stage of magma ocean evolution. An 
ilmenite-bearing cumulate (IBC) layer would begin 
to form after about 95% solidification of the magma 
ocean [23]. This means that the magma ocean cumu- 
lates were stratified with dense FeO-rich phases 
overlying less dense MgO-rich phases. 

This configuration of more dense above less dense 
mantle is gravitationally unstable. This paper ex- 
plores the chemical and thermal consequences of an 
internal evolution model accounting for the possible 


role of this chemical buoyancy. Large scale overturn 
of the lunar interior would redistribute the compo- 
nents of the mantle with interesting petrological 
[2,3,15,24,25] and thermal consequences. Since it 
formed from the last remaining liquid in the magma 
ocean, the IBC has a high content of incompatible 
radioactive elements. After overturn, radioactive 
heating within the IBC heats the overlying mantle, 
causing it to melt. This heating results in deep 
melting and, for some parameters, melting occurs 
after a reasonable time interval to explain the onset 
of mare basalt volcanism. If the Moon formed by a 
giant impact into the Barth [26] and consists primar- 
ily of material derived from the Earth's mantle, 
material for a metallic-Fe core may not be present. If 
there is no metallic core, IBC would sink to the 
center of the Moon thus providing a simple way to 
create a relatively dense core. If a metallic core were 
already present, the IBC would form a layer sur- 
rounding the core. 

2. Solidification of the magma ocean: initial chem- 
ical stratification 

We assume that a magma ocean differentiated 
into an anorthositic crust and a chemically stratified 
mantle. The products of this differentiation depend 
on the depth of the magma ocean, the lunar bulk 
composition, and the nature and efficiency of the 
differentiation process. Estimates of magma ocean 
depths range from 250 km to the entire Moon [27]. It 
has been argued that the absence of global scale 
extensional or compressions! features on the lunar 
surface limit variations in the lunar radius to ± 1 km 
over the past 3.8 Gyr [9]. Solomon and Chaiken [9] 
found that for this radius constraint to be satisfied, 
the contraction due to cooling of the outermost lay- 
ers of the Moon, which they assumed to be the 
crystallized product of the magma ocean, must be 
balanced by radioactive heating of the deeper inte- 
rior. They argued that this constraint limits magma 
ocean depths to less than 400 km. Kirk and Steven- 
son [12] indicated that a 630 km deep magma ocean 
can be allowed if contractional cooling of the magma 
ocean is partially offset by volume changes due to 
partial melting and chemical differentiation. 

These magma ocean thickness contraints pose an 
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obstacle to theories of lunar origin that call lor an 
initially molten Moon [10], But other taclors not 
considered in previous models may be important. 
The downward transport of heat-producing elements, 
as discussed later, would hasten the cooling ol the 
shallow Moon in the period preceding 3.K (ia and 
increase the heating of the deep interior on longer 
time scales. Both of these effects would allow a 
larger magma ocean thickness to be consistent with 
the Solomon and Chaiken constraint. 

A magma ocean shallower than 400 km is prob- 
lematical from a chemical perspective. Unless the 
bulk lunar A1 2 0 3 content is much higher than about 
4% [28-30], a magma ocean thinner than 400 km 
cannot provide the A1 2 0 3 content of the anorthositic 
crust [18]. Although we cannot rule out an AJ 2 0 3 -en- 
riched Moon, we adopt the generally assumed bulk 
Moon composition. 

Recognizing these major uncertainites in deter- 
mining the thickness of the magma ocean, we con- 
sider values ranging from 400 to 800 km. Our stan- 
dard case is an 800 km deep magma ocean that 
creates 60 km of anorthositic crust (including later 
additions of Mg* suite plutons), 30 km of an il- 
menite-enriched cumulate layer beneath the crust, a 
710 km cumulate layer with relatively dense (low 
Mg*) mafic orthopyroxene-olivine layers near the 
top and less dense (high Mg * ) olivine-rich cumulate 
layers near the base (e.g., [31]). This stratigraphic 
sequence overlies 940 km of primitive Moon (see 
[30] for a discussion of chemical models of the 
Moon). 

The ilmenite liquidus in mafic to ultramafic sili- 
cate melts occurs for melts containing about 4-5 
wt% TiO z at temperatures less than 1100°C and at 
low pressures [23]. If we assume a Ti0 2 content for 
the magma ocean of 0.2-0.3 wt% Ti0 2 [30] then 
IBC forms a layer about 30 km thick. The density of 
the IBC layer is controlled by the content of ilmenite 
( p ~ 4800 kg/m 3 ), the Mg’ values of pyroxene 
and the modal abundance of plagioclase ( p ~ 2750 
kg/m 3 ). The normative ilmenite content of ilmenite 
saturated melts at about 1100°C is about 10 wt% 
[23,32] and the normative content of plagioclase is 
about 35 wt%. If a significant fraction of the plagio- 
clase floated and was added to the anorthosite crust, 
then the ilmenite content of the cumulate layer could 
be as high as 15 wt%, and the thickness of this layer 


would be about 20 km. We assume very little modal 
plagioclase in t he cumulate layer and a range ot 
reasonable Mg* values for the modal pyroxene 
(Mg ’ = 0.2-0.5). This range of Mg * values is con- 
sistent with the Mg* values of the most evolved 
ferroan anorthosites (Mg* ~ 0.40) which have crys- 
tallized from silicate melts (presumably the magma 
ocean) that have not reached the ilmenite liquidus 
[27,33]. The density of a cumulate with 15% ilmenite 
and 85% pyroxene (Mg* =0.30) is about 3700 
kg/m 3 

The most primitive mafic minerals in ferroan 
anorthosite have Mg * = 0.76 [24]. If the ferroan 
anorthosite appeared after 80% crystallization of the 
magma ocean [22] then the cumulates lying below 
the IBC layer must include cumulates with Mg * > 
0.76. These cumulates are likely to be largely 
harzburgites (see [31] for review) with p ~ 3400 
kg/m 3 . The first cumulates at the base of the magma 
ocean cumulate layer are likely to be dunites with 
Mg* values determined by the compositon of the 
magma ocean and the pressure of crystallization [34]. 
The Mg* values of olivine in equlibrium with a 
peridotitic magma with Mg * > 0.80 at pressures of 
3-4 GPa are Mg * > 0.90 [34]. These cumulates 
have p-3300 kg/m 3 . A calculated initial density 
distribution for the cumulate layer consistent with 
these chemical constraints is given in Fig. 1. 

The time needed to solidify the magma ocean is 



Fig. 1. Intrinsic (chemical) density distribution as a function of 
radius before (“undifferentiated ,t ) and after (‘‘differentiated ) 
overturn. We assume that buoyancy forces in the original unstable 
configuration result in a rearrangement of mantle layers so that 
density decreases monotonically with radius. 
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largely a function of the thickness and persistence of 
a solid layer at the surface through which heat is 
transferred only by conduction. Prior to the forma- 
tion of a stable solid crust, convective cooling and 
solidification of the magma ocean will be rapid (cf. 
[19,20]). The rate of solidification would slow dra- 
matically, however, once the magma ocean is blan- 
keted by a stable plagioclase-rich crust [16]. The rate 
of solidification should be about 100 km per 10 8 yr 
after a stable solid crust of 10 km or more is a 
permanent feature of the magma ocean. According to 
Phinney [22] the anorthosite crust reached this stage 
after approximately 80% of the magma ocean had 
solidified. 

Thus, the earliest cumulates, precipitated from a 
magma ocean undergoing rapid convective cooling, 
formed within a short time interval of a few million 
years. The thickness of this cumulate sequence is 
approximately 700 km (assuming the anorthosite 
crust forms only after about 80% solidification of the 
magma ocean) for an 800 km thick magma ocean 
(e.g., [35]). The last 100 km of the magma ocean, 
crystallized over an interval of roughly 100 Myr, 
fractionated into about 40 km of anorthosite crust 
(assuming 20 km of the present-day crust is com- 
posed of younger Mg * suite plutonic rocks [36]) and 
about 60 km of FeO-enriched pyroxenite cumulates, 
the upper 20 km of which are the IBC layer. As 
discussed later, the time needed to form the IBC 
layer and the Iate-stage pyroxenite cumulates ex- 
ceeds by one to two orders of magnitude the time 
scales during which the lowermost olivine-rich cu- 
mulate sequence could establish a gravitationally 
stable stratification. Therefore, the IBC layer formed 
above a relatively stable platform of MgO-rich, dunite 
cumulates, which themselves overlay olivine -ortho- 
pyroxene cumulates, representing an inverted cumu- 
late stratigraphy. 

3. Large scale overturn of initial unstable stratifi- 
cation 

The cumulate mantle resulting from solidification 
of the magma ocean is gravitationally unstable with 
dense ilmenite-rich cumulates overlying olivine-or- 
thopyroxene cumulates with Fe/Mg that decreases 
with depth. This density stratification will cause 


large scale overturn of the lunar interior. C onsider 
first a simple idealization that ignores mixing during 
overturn. Primitive lunar mantle rises to its level of 
neutral buoyancy in the cumulate layer and most of 
the dense ilmenite-rich layer sinks to the center, 
forming a core [15]. The mean density of the IBC 
(p~ 3700 kg/m 1 with 15% ilmenite and 85% 
clinopyroxene with Mg* =0.3) is significantly 
greater than the density of the lower cumulates and 
the undifferentiated Moon, assumed to be garnet- 
Iherzolite ( p ~ 3350 kg/m 3 ). For an 800 km thick 
magma ocean, the IBC, after sinking through the 
mafic cumulates and undifferentiated interior, segre- 
gated at the core-mantle boundary as a 200 km thick 
layer if a 400 km radius Fe-Ni-S core already 
existed or as a sphere of about 550 km radius, if 
there was no pre-existing core. An idealized stratig- 
raphy and density distribution after the establishment 
of a gravitationally stable lunar mantle is shown in 
Fig. 1. The moment of inertia of the Moon, which 
only slightly exceeds that for a uniform density 
sphere (cf.[37]), can be explained by a small ( ~ 400 
km radius) metallic core or by a model with no 
metallic core in which the mantle density increases 
with depth as implied by the model in Fig. 1 . A 400 
km thick magma ocean, rather than 800 km as 
assumed in this example, would result in an IBC 
core radius of about 475 km. 

After the overturn the uppermost mantle would 
probably still contain a thin zone perhaps less than a 
few kilometers in thickness of incompatible 
element-rich cumulates crystallized from the last liq- 
uids of the magma ocean (which have been termed 
UrKREEP cumulates; cf. [38]) and whatever fraction 
of the IBC layer that survived the overturn. These 
cumulates rest upon olivine-cumulates (Mg * > 0.90) 
which, in turn, overlie a 150 km thick layer derived 
from the initially undifferentiated interior. This layer 
was displaced upward by the sinking IBC. The layer 
corresponding to the undifferentiated interior over- 
lies a cumulate section which is initially olivine-rich 
but grades with depth to an olivine-orthopyroxene 
assemblage. The olivine-pyroxene assemblage con- 
tinues almost to the IBC layer-mantle boundary. 

The true stratigraphy of the Moon is bound to be 
more complex since at least some mixing will occur. 
Small masses of IBC, for example, have small set- 
tling velocities and may remain trapped and dis- 
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persed throughout the mantle. Sinking and rising 
diapirs may entrain and conveetively mix with por- 
tions of the mantle through which they pass [3,39]. 
The sinking and associated mixing of cumulates 
represents a process by which KREEP-rich cumu- 
lates and IBC could be dispersed throughout the 
mantle (see also [19,40]). Such mixing would reduce 
the density stratification shown in Fig. I, which 
assumes no mixing at ail. 

The roughly 550 km radius IBC core contains 
about 15 wt% ilmenite. The IBC cumulates should 
also be enriched in highly incompatible elements, 
especially the heat-producing elements. The IBC 
were derived from melts which contained about 20 
times lunar abundances of heat-producing elements 
because the melts parental to IBC are liquids ob- 
tained after 95% crystallization of the magma ocean. 
We assume that rising plagioclase cumulates are 
impoverished in these elements, so that the comple- 
mentary mantle contained about 30 times lunar abun- 
dances. Some fraction of these heat-producing ele- 
ments will be carried to the lunar IBC core. We 
assume that the heat-producing elements range from 
5 to 10 times bulk lunar abundances in the analysis 
that follows. 


3.1. Temporal and spatial scales for overturn of 
initial chemical stratification 


Rayleigh-Taylor instability will take two forms, 
large scale overturn of the interior driven by the high 
density of the ilmenite-rich cumulate layer and inter- 
nal overturn of the magma ocean cumulates beneath 
the ilmenite-rich cumulate layer. The magma ocean 
cumulates that form below the IBC layer are inter- 
nally stratified with denser Fe-rich cumulates overly- 
ing earlier formed Mg-rich cumulates. Most of these 
cumulates are less dense than the primitive interior 
and so do not promote large scale overturn of the 
lunar interior. The time scale for internal overturn of 
the cumulates can be estimated from the growth rate 
of Rayleigh-Taylor instability in a constant viscosity 
fluid layer in which the fluid density varies linearly 
with depth: 


A = 2 d 

where d is (he depth of the layer, dp is the density 
difference between the top and bottom, and p is the 
viscosity. For dp - 350 kg/m\ p = l() n Pa s and 
</ = 7X10* m, the characteristic time is roughly 
0.03 Myr. We argued previously that the crystalliza- 
tion of the last 100 km of the magma ocean would 
slow dramatically once an anorthosite crust was sta- 
bilized. Estimates of the crystallization time are on 
the order of 100 Mvr. It appears, therefore, that the 
lower cumulate pile will have reached gravitational 
equilibrium before most of the IBC has formed. 

Large scale gravitational instabi It y will be driven 
primarily by the high density of the IBC layer. First, 
consider the simple, but probably unrealistic, case in 
which the IBC solidifies rapidly compared to the 
time scale of differentiation. The time scale for the 
sinking of the IBC layer to the center of the Moon 
can be controlled by the time for diapirs to form by 
Rayleigh-Taylor instability or by the time for the 
resulting diapirs to sink through the lunar mantle. 
This time scale will be determined by the wave- 
length of instability and the corresponding size of the 
diapirs. 

An estimate of the time needed for Rayleigh- 
Taylor instablity of the ilmenite-rich cumulate layer 
to develop is (e.g., [52]): 


6.5m 2 /V/ 


Apgh 

with a wavelength of instability: 


A = 2.9 h 


where p x /x 2 . In these expressions p { and p 2 are 
the viscosity of the IBC and underlying cumulates, 
respectively (1G 19 Pa s for mantle silicates near their 
melting temperature would be a reasonable minimum 
value), Ap (0.35 X 10 3 kg/m 3 ) is the density differ- 
ence between them, h is the layer thickness ( ~ 20 
km), and g is the acceleration of gravity at the 
appropriate depth. Given the uncertainty in a number 
of these parameters, particularly the viscosities, this 
time ranges from 10 4 to 10 5 yr. The wavelength of 
instability is sensitive to the viscosity of the il- 
menite-rich cumulate layer. If this viscosity were an 


r 
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order of magnitude less than that of underlying 
cumulates, ~ 10 t,lc wavelength of insta- 

bility would be about 6 h or 120 km. Equating the 
layer volume in one wavelength of instability, A ~/j, 
to that of the resulting diapir gives a diapir diameter 
of about 80 km. About 2000 such diapirs would 
form. Such a small scale of instability could result in 
some mixing during the differentiation. This number 
of diapirs varies inversely with A 2 and therefore the 
number of diapirs is proportional to ( /x,/ /x 2 ) 2/ \ For 
/x ,// lx 2 = 10” 5 , only a single diapir would form. 
This larger scale of instability would result in less 
mixing during the differentiation. It is interesting to 
hypothesize that the lunar center of figure -center of 
mass offset and corresponding crustal asymmetry (cf. 
[37]) might be indicative of such a large scale insta- 
bility. 

The time required for a diapir to sink is estimated 
using the Stokes settling velocity of a fluid sphere 
with viscosity /lx, settling in a fluid of viscosity yx 2 
(cf. [41]): 

v _ 1 | 
e -2 1 1 1 ! 2^1 / 

where a is the diapir radius 4tt(i 3 /3 = \ 2 h. The 
sinking time over a distance comparable to the radius 
of the Moon ( ~ 1740 km; see Fig. 1) for the cases 
with /x,//z 2 ~ 10” 1 and 10 s cited above with /x 2 
= 10 19 Pa s are 2 Myr and 3 X 10” 2 Myr, respec- 
tively. Thus for jxj/^ « 10” 1 the differentiation 
time is controlled by the sinking of small diapirs but 
the opposite is true for /x,/^x 2 « 10” 5 . Good esti- 
mates of the relative viscosities are thus important 
for understanding the rate controlling step in the 
differentiation. Whether the IBC layer is much less 
viscous or only slightly less viscous than the under- 
lying mantle depends strongly on the melt fraction 
present as diapiKs) begin to form and on how this 
melt is distributed on the grain-size scale. If melt 
wets grain faces forming a thin, low shear strength 
film, a small melt fraction may greatly reduce the 
viscosity. Alternatively, if melt does not wet grain 
faces, small amounts of melt along grain edges and 
at grain triple junctions will only modestly reduce 
the viscosity (cf. [42D. Experimental evidence on the 
degree to which melt of the appropriate composition 
wets grain boundaries of coexisting minerals is not 
yet available. 


it is important to recognize that the differentiation 
process will be much more complex than these sim- 
ple estimates appear to imply. Diapirs will develop 
while the layer is still solidifying, rather than after 
solidification is complete as assumed above. The 
scale of diapirs that develop can be estimated by 
equating the thickening rate of the solidified IBC 
layer k/h with the rate of diapir formation as deter- 
mined by the reciprocal of the Rayleigh-Taylor 
growth time given above. This gives an estimate of 
the steady-state layer thickness as solidification pro- 
ceeds: 

6.5p}/ 3 p\ /3 h \ 

Apg ) 

where Ap is density difference between the layer 
and underling cumulates and other symbols have the 
same meaning as above. For this layer thickness an 
estimate of the wavelength of instability and the 
resulting diapir size also follow from expressions 
given above. For /tx,//u, 2 ~ 10 1 and a solidification 
rate of 0.1 km/Myr, the estimated layer thickness, 
wavelength of instability, and diapir diameter are 0.5 
km, 3 km, and 2 km, respectively. For /Xj//x 2 — 10” 5 
and the same solidification rate, the corresponding 
values are 0.1 km, 13 km, and 4 km. 

The rate of solidification of the layer appears to 
be an important parameter. The Solomon and Longhi 
[14] model predicts that the residual layer remains 
partly liquid for 350 Myr after crystallization of the 
majority of the magma ocean. This gives a solidifica- 
tion rate of about 0.05 km/Myr. In their model solid 
mantle is present above the residual liquid layer. 
Sinking of this gravitationally unstable layer through 
the residual layer approximately halves the distance 
through which heat must be conducted to reach the 
surface (see [14], fig. 8). The solidification time of 
the layer by steady state heat conduction is propor- 
tional to the layer depth (the distance over which 
heat must be conducted) and so would be reduced by 
a factor of about 2. For transient heat conduction the 
cooling time is proportional to the square of the 
depth and so would be reduced by a factor of 4. 

Diapir size is proportional to {k so that even with 
solidification rate 2-4 times larger, the diapir sizes 
will be small. 

Thus, as the IBC layer crystallizes, small diapirs 
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will continually form and sink into the underlying 
cumulates generating a layer that is a rmxlure of IBC' 
and earlier cumulates. The thickness of this mixed 
layer can be estimated as the product of the solidifi- 
cation time and the diapir settling velocity. For a 20 
km thick IBC' layer solidifying at 0.1 km/Myr (so- 
lidification time = 2(H) Myr), a 2 km diameter diapir 
should sink a maximum distance ot 80 km. The IBC 


layer is thus diluted by a factor of 4. Since the 
settling velocity is proportional to h and the solidifi- 
cation time is inversely proportional to h , the mixed 
layer thickness is actually independent of the solidi- 
fication rate. The size of IBC diapirs in this mixed 
layer is, however, dependent on this rate. 

This thicker mixed layer will generate global 
scale instability as initially described above for the 
IBC layer alone. The time scale of instability, which 
depends only the product of Ap and layer thickness, 
remains the same as that estimated previously. The 
wavelengths of instability will be larger by about a 
factor of 4 corresponding to the greater thickness of 
the mixed layer. However, this mixed layer itself 
will form diapirs on time scales shorter than the time 
for the whole layer to form. The differentation pro- 
cess is clearly complex, but these simple considera- 
tions suggest that differentiation will occur quickly 
and a dilution of IBC by about a factor of 4. 

How heat sources will be distributed in this sink- 
ing IBC-cumulate mixture is more difficult to assess. 
If heat-producing elements were perfectly incompati- 
bile, they would be continually concentrated into 
whatever liquid remained, until that liquid became 
saturated in these elements. Sinking solid will, how- 
ever, carry with it some buoyant intragranular liquid. 
Adiabatic compression will freeze this liquid increas- 
ing its density and allowing it to continue sinking. 
How much liquid can be carried downward in this 
way depends on the rate of melt migration and/or 
the rate of compaction of the solid. Given the obvi- 
ous complexities of the differentiation process, we 
simply consider a range of heat-production values in 
the IBC core from 5 to 10 times bulk Moon values. 

3.2. Internal evolution after large scale overturn 


The initial chemical density stratification of the 
lunar mantle after overturn is large and gravitionally 
stable (Ap~ 200 kg/m 3 ; see Fig. 1). For a typical 


coefficient of thermal expansion (3 X 10 * /K), the 
supcradiabatic temperature difference between the 
top and bottom of the lunar mantle would have to 
exceed 2(KH)°C for large scale thermally-driven man- 
tle convection to disrupt this initial stratification, 
liven assuming very low initial temperature in the 
primitive lunar interior beneath the magma ocean, 
known crystallization temperatures of the magma 
ocean would not allow an initial temperature differ- 
ence this large. Thermally-driven convection in the 
interior will begin in mixed layers near the top of the 
mantle and at the core-mantle boundary with thick- 
nesses determined by the chemical density difference 
that can be penetrated by thermally buoyant mantle. 
The following analysis treats the development of 
such a layer adjacent to the core-mantle boundary 
assuming that the effect of the complementary layer 
expected to develop at the top of the mantle has not 
yet reached a sufficient depth to be important. 

With a concentration of incompatible heat-produc- 
ing elements 5—10 times that of the bulk Moon, the 
ilmenite core heats the overlying olivine-orthopy- 
roxene cumulates. A thermal boundary layer at the 
core-mantle boundary that thickens by conduction 
becomes unstable when the Rayleigh number based 
on the boundary layer thickness exceeds a critical 
value: 

cnl 

ILK 

generating plumes that rise into the overlying chemi- 
cally layered mantle. Here T c and 7^ are tempera- 
tures at the core mantle boundary and in the overly- 
ing mantle, a is the coefficient of thermal expansion 
(3 X 10” 5 K -1 ), k is the thermal diffusivity (10“ 6 
m 2 /s), and g is the gravity at a depth of 1200 km. 
This instability criterion assumes that the chemical 
density difference over the thickness of the thermal 
boundary layer is negligible compared to the thermal 
density difference. The time r for the thermal 
boundary layer to grow conductively to thickness 8 
at which it becomes unstable is: 

_ 8 2 

2k 

Repeated growth and instability of the boundary 
layer forms thermals that rise into the overlying 
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mantle, faking (7^ + T m )/2 as an estimate of aver- 
age temperature of material in the thermal boundary 
layer, the heat flux ( q ) into the overlying mantle is 
simply the product of the amount of heat added each 
time the boundary layer becomes unstable pc (T c - 
7 m )5/2 and the frequency of instability 1/r: 

P^(T C -T m ) 8 
<7 = o — 

2 7 

where c p (1 kJ/kg K) is the heat capacity. 

Several simple but reasonable assumptions are 
made to treat the dynamics of thermal plumes rising 
in chemically stratified surroundings. (1) Plumes 



thermal expansion: 



Density 


neutral buoyancy level 



Fig. 2. Important features of the thermal evolution model of a 
chemically stratified lunar mantle heated from below by radioac- 
tivity in the ilmenite core. The bottom diagram illustrates that 
plumes form from an unstable thermal boundary layer at the base 
of the mantle. Plumes rise adiabatically to a height at which the 
combined themai and chemical buoyancy in the plume and in 
overlying mantle are equal. The rise of plumes is assumed to 
create a thermally and chemically well-mixed layer. At any time t, 
plumes, which rise through the well-mixed layer of thickness 
**(D, penetrate the overlying mantle to a height z'(f+r), 
where r is the time required to conductively form an unstable 
thermal boundary layer. The mantle is taken to have an initial 
temperature Tj with an intrinsic (chemical) density that decreases 
linearly with height, as estimated from Fig. 1. 



Fig. 3. Calculated thermal structure as a function of time for a 
core radioactivity 10 times bulk earth, a chemical density stratifi- 
cation above the core-mantle boundary of 0.9 X 10~ 7 /km, a 10 19 
Pa s mantle viscosity, and an initial core and mantle temperature 
of 1200°C. This corresponds to parameters in the third row of 
Table 1 and the results labelled (3) in Fig. 4. 


formed at the core-mantle boundary (CMB) rise into 
intrinsically (chemically) less dense but colder man- 
tle. Plumes with the average temperature and compo- 
sition in the marginally stable thermal boundary 
layer rise to a height above the CMB where the 
positive thermal and negative chemical buoyancy 
balance. (2) Plumes cause overturn and mixing of the 
mantle through which they rise, thus thermally and 
chemically homogenizing a mantle layer equal to the 
plume rise height. Thus, as shown in Fig. 2, the 
mixed layer above the thermal boundary layer has a 
uniform temperature T m . The initial intrinsic chemi- 
cal density stratification is taken to be a simple linear 
function of height above the CMB. With the assump- 
tion of complete mixing, the intrinsic density in the 
mixed layer is just the average of the initial, pre-ex- 
isting variation over the layer thickness, as shown in 
Fig. 2. Hie density of plumes and of the mantle 
above the mixed layer are given by the intrinsic 
density reduced by thermal expansion, as shown in 
Fig. 2. As the temperature T m increases with time, 
the mixed layer thickness also increases. 

The evolution of the mixed layer thickness z* 
and the mean temperature T m within it follow from a 
simple energy balance: the change of thermal energy 
over a small increment of time is a consequence of 
heat added by thermal boundary layer instability and 
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entrainment of cooler overlying mantle into the mixed 
layer. In a spherical geometry, where R is the radius 
of the CMB, this energy balance can be written: 

dz’ 

= q -pc p (T m -T t ) — 

The core temperature is similarly determined by 
balancing secular and radioactive heating with the 
rate of heat loss to the overlying mantle. Since 
substantial melting of the ilmenite-rich core and the 
thermal boundary layer occurs, the heat of melting 
(600 kJ /kg K) is also included. 

We explore predictions of this model for a range 
of mantle viscosity (10 19 -10 :o Pa s) and core ra- 
dioactivity values (5-10 times bulk earth) with an 
initial core and mantle temperature of 1200°C. This 
initial temperature of the core and the mantle just 
above it is estimated by taking the solidification 
temperature of the IBC layer to be 1100°C and by 
including the increase in temperature accompanying 
the adiabatic compression of the sinking IBC 
diapiKs). Just above the IBC core in Fig. 1, the 
density changes by about 1% over a height of 100 
km. However, if some mixing occured during over- 
turn, the compositional gradient would be reduced. 
We therefore consider cases with 0.9%, 0.6%, and 
0.3%, density changes over 100 km. Table 1 summa- 
rizes the range of model parameters the we have 



Fig. 4. Mixed layer thickness and temperature as functions of time 
for a range of chemical density gradients corresponding to rows 
(1), (2X and (3) in Table 1. The case with the highest chemical 
density gradient (3) corresponds at a given time to the thinnest 
mixed layer and the highest temperature, since nearly the same 
amount of heat has been added to a thinner layer. 


Table 1 

Selected model parameters and results 


heating 

viscosity 

(Pa-sec) 

cote 

radius 

(km) 

compositional density 
gradient 
(lcm->) 

time to 
melting 
(Myr) 

melting 

depth 

(km) 

melting 

temperature 

(°Q 

melting 

pressure 

(GPa) 

10. 

10 19 

550. 

0.3xl0- 7 

996. 

336. 

1365. 

1.84 

10. 

10 19 

550. 

0.6xl0- 7 

636. 

760. 

1458. 

3.08 

10. 

10> 9 

550. 

0.9x10-2 

507. 

807. 

1491. 

3.54 

5. 

10' 9 

550. 

0.3x10-2 

1965. 

522. 

1421. 

2.59 

10. 

1020 

550. 

0.9x10-2 

735. 

664. 

1459. 

3.09 

10. 

10' 9 

475. 

0.6x10-2 

774. 

810. 

1474. 

3.42 

10. 

10 19 

475. 

0.9x10-2 

642. 

948. 

1502. 

3.79 
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examined. The pendotite solidus is taken to he that 
of an ilmenite-bearing har/burgitc [3]. Fig. 3 shows 
calculated temperature distributions lor one of the 
models. The temperature distribution within the ther- 
mal boundary layer, approximated as a simple linear 
variation, is shown just prior to instability at each of 
the times shown. Fig. 4 shows the variation in mixed 
layer temperature and thickness with time for several 
of the models listed in Table I. These are further 
discussed below. 


4. Discussion 

In all of these models melting in the mantle 
begins within the thermal boundary layer just above 
the core-mantle boundary. Consider, for example, 
the thermal evolution of the mantle for model (3) 
(Table 1) with a heat source of 10 times bulk earth 
values and a density gradient of —0.9% per 100 km 
in height above the core-mantle boundary. The tem- 
perature profile of the mantle above the core-mantle 
boundary 400 Myr after the overturn of the lunar 
mantle and core formation consists of a 320 km 
well-mixed mantle layer with a temperature of 
1425°C, a temperature still well below the mantle 
solidus at that depth (Fig. 3). However, the roughly 
70 km thick, thermal boundary layer is partially 
molten. Melting in the thermal boundary layer may 
slightly increase or decrease the density of the par- 
tially molten mantle, given estimates of melt and 
solid densities at these conditions [5]. The effect of 
melting on buoyancy is likely to be small compared 
to the effect of thermal expansion. We therefore have 
neglected melt buoyancy in the models described 
above. The possible effect of the presence of melt in 
the lunar core will be considered later. 

The thermal boundary layer is thermally buoyant 
and becomes unstable, mixing both thermally and 
chemically with the overlying mantle. The mass of 
the thermal boundary is cooled to the mean tempera- 
ture of the mixed mantle. With the temperature of 
the mixed layer below the local solidus, whatever 
melt is transported to the mixed layer from the 
thermal boundary layer crystallizes and releases la- 
tent heat in the mixed mantle layer. The temperatures 
depicted in Fig. 3 include the thermal contribution of 
this latent heat. 


Within 500 Myr of core formation, the thickness 
of the mixed layer has reached over 340 km, and 
more significantly, the temperature at the top of the 
mixed layer has exceeded the local solidus at a depth 
corresponding to a pressure of 3.5 GPa. The timing 
and depth at which melting begins in the mixed layer 
is dependent on several of the model parameters. 
Fig. 4a and b describes the time evolution of mixed- 
layer temperature and thickness, respectively, for 
density gradients of —0.3, —0.6 and —0.9%/ 1 00 
km in height above the core-mantle boundary. The 
coolest and thickest mixed layer occurs for the man- 
tle with the smallest density gradient. Greater thick- 
ness for smaller density stratification is a conse- 
quence of the greater rise height of thermal plumes 
in a smaller adverse compositional density gradient. 
Lower temperature is a consequence of adding heat 
from repeated instabiltiy of the thermal boundary 
layer to the larger mass of mantle in the mixed layer. 
Based on the results of Fig. 4, the mixed mantle with 
a low compositional density gradient at a given time 
is thicker but cooler than that with a high composi- 
tional density gradient. The stable chemical stratifi- 
cation acts to confine heat to a thinner layer. As a 
consequence, melting in the mixed layer occurs at 
lower pressure, lower temperature and more than 400 
Myr later for the low density gradient than in the 
high density gradient (Table 1). 

The thermal evolution and subsequent melting of 
the mixed layer also depends on the concentration of 
heat-producing elements in the core. Reducing the 
radioactivity to 5 times bulk earth values has dra- 
matic consequences. The solidus at the top of the 
mixed layer is reached but only after a considerable 
time interval 1 Gyr later than the case for 10 times 
core radioactivities (Table 1). Similarly, increasing 
the mantle viscosity from 10 19 to 10 20 Pa s increases 
the time for melting at the top of the mixed layer by 
about 200 Myr. Melt is generated at slightly lower 
temperatures and somewhat shallower (compare (3) 
and (5) in Table 1). 

We already noted that the temperatures obtained 
at the bottom of the thermal boundary and in the 
core exceed the solidus soon after the overturn. 
Assuming that the mantle is a harzburgite or perhaps 
an olivine-bearing orthopyroxenite as seems likely, 
the melt fraction can approach 100%. The tempera- 
ture of the forsterite-enstatite (Mg* *= 1.0) cotectic 


/> ( i/ CSSt i; \f l\irmcnlH'r / i.iuth and Waneuiry Science Letters l -U l IW5) >01 


5 1 I 


ckness 
ii, and 
of the 
^ depth 
liming 
:d layer 
meters, 
mixed- 
;ly, for 
%/l00 
ry. The 
ic man- 
r thick- 
conse- 
plumes 
;radient. 
mg heat 
oundary 
;d layer, 
tie with 
en time 
omposi- 
stratifi- 
er. As a 
ccurs at 
than 400 
n in the 

siting of 
Oration of 
icing the 
has dra- 
>p of the 
siderable 
10 times 
ncreasing 
increases 
l layer by 
ttly lower 
npare (3) 

, obtained 
nd in the 
overturn. 
>r perhaps 
ms likely, 
i tempera- 
)) cotefctic 


is about 19()()°C at 4 GPa; the congruent melting 
points for enstatite and lorsterite are about 1920 C 
and 2100°C at these pressures [43]. Hie Mg * values 
of these cumulates is assumed to be 0.80. I lie liq- 
uidus temperature ot olivine is depressed by 
3°C/ Mg ‘-number or by about 60°C for an Mg* = 
0.80 (see [44], fig. 16.2). Thus, the liquidus tempera- 
tures of olivine and enstatite (Mg * = 0.80) is roughly 
2040°C and 1860°C, respectively, at these pressures. 
Thus for the temperatures shown in Fig. 3, the lower 
portions of the thermal boundary layer are totally 
molten for pyroxene-rich lithologies but only par- 
tially molten for olivine-rich cumulates. 

The liquidus temperature for the lunar core is 
even lower because the Mg* values are lower. If 
this analysis is correct, the lunar core is totally 
molten about 400-500 Myr after core formation for 
10 times bulk earth radioactivities. It is interesting to 
note that the high apparent Love number of the 
Moon can be explained by a partially molten central 
region with 800 km radius that persists to the present 
[13]. The concentration of heat sources and low 
melting temperature materials toward the center of 
the Moon would make it easier to explain a currently 
molten lunar core. 

What role might melting of the lunar IBC core 
play in the thermal and chemical evolution of the 
lower mantle? In our model, the lunar core contains 
about 15 wt% ilmenite and 85 wt% pyroxene ± 
plagioclase. The density of the molten core can only 
be estimated given the large extrapolations needed to 
obtain the appropriate thermal expansions and com- 
pressibilities. Moreover, reasonably modest changes 
in the Ti0 2 , Mg* number and the Si0 2 content can 
have significant effects on the calculated density. For 
mare-like melts with 8% Ti0 2 (equivalent to 15% 
normative ilmentite), the results of Delano [5] sug- 
gest that the density of the molten core is about 
3400 ± 100 kg/m 3 . 

In comparison, the density of a partially or even 
fully molten mantle in the thermal boundary layer is 
much less. We assume that the mantle in the thermal 
boundary layer is largely low CaO-pyroxene (Mg 
ss 0.80) ± olivine and minor ilmenite. The density of 
molten low CaO-pyroxene at 1900°C is about 2700 
kg/m 3 [45], whereas any unmelted olivine has a 
density of less than 3300 kg/m 3 at these tempera- 
tures [5]. The density of the melt corrected to 4.5 


GPa [5] is roughly 3000 kg/m 3 . It follows that the 
partially to totally molten thermal boundary layer is 
gravitationally stable relative to the molten core. The 
two layers may not remain isolated from each other, 
however. Core material may be entrained as the 
thermal boundary layer becomes unstable and rises 
in the overlying cumulates; alternatively convection 
in the core may entrain melted mantle in the thermal 
boundary layer, eroding the base of the mantle. 
These processes may be important but have not been 
incorporated in the model described above and are 
not considered further in this paper. Both effects 
would enhance the rate of heat transfer into the base 
of the mantle from the core. 

If our lower estimates for the densities of the 
molten core are correct, it may be unstable relatively 
to the overlying mantle. A molten core with p < 3400 
kg/m 3 is capable of becoming unstable relatively to 
the more iron-rich pyroxene-rich cumulates (Fig. 1). 
If so, core material, once it melts, might rise to form 
a layer displaced from the center of the Moon, at its 
height of neutral buoyancy in the overlying cumulate 
mantle. It also is possible that the core material 
might mix with the overlying mantle to a limited 
extent. These possibilities, although interesting, are 
not pursued further in this paper. 

4.L Mare basalt genesis 

An important constraint on mare petrogenesis pro- 
vided by experimental phase relations shows that 
pristine mare glasses are multiply saturated with 
olivine and orthopyroxene at pressures between 2.0- 
2.5 GPa but have only olivine on the liquidus at 
lower pressures [4,46]. If they are primary magmas 
and if both olivine and orthopyroxene were left in 
the source, the mare glasses must have been derived 
from pressures of about 2.0-2.5 GPa. The pressures 
would be lower and higher than this range if only 
olivine or orthopyroxene, respectively, were left in 
the source. We assume that the picrite mare glasses 
are primary magmas derived by relatively low de- 
grees of melting [47]. The source rocks in our model 
were derived as olivine-orthopyroxene cumulates 
from a shallow magma ocean and then carried down- 
ward to higher pressure in the interior during the 
overturn. For pressures greater than about 0.5 -0.6 
GPa, the olivine-orthopyroxene cotectic lies very 
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close to the plagioclase— orthopyroxene join and the 
cumulates deposited from these liquids would have 
orthopyroxene olivine. It follows, therefore, that 
olivine cannot be the sole unmeited phase in the 
source, which means that the mare glasses were 
derived from an orthopyroxene ± olivine source, and 
must have been derived from depths corresponding 
to pressures of at least 2.0— 2.5 GPa ( ~ 400—500 
km). 

If the mare basalts are not primary magmas but 
have undergone near-surface crystallization of 
olivine, then the estimated depths of generation are 
greater than the pressures of multi-saturation. This 
follows from the experimental observation that the 
orthopyroxene-olivine cotectic becomes more 
olivine-normative with increasing pressure. If the 
mare basalts are collections of polybaric melts as 
proposed by Longhi [481, then the apparent pressures 
of multi-saturation represent a mean pressure of gen- 
eration, i.e., melts were generated at pressures both 
above and below 2—2.5 GPa. 

As a source for mare volcanism, melting at the 
top of the mixed layer must occur at low enough 
pressure for melt to be buoyant relative to solid and 
at high enough pressure to satisfy the multi-satura- 
tion contraint discussed above. Melting at the top of 
the mixed layer must begin within 500-600 Myr of 
the overturn since the oldest high-Ti0 2 mare basalts 
appeared at lea^t 3.9 Myr ago [6]. Our results sum- 
marized in Table 1 show that core radioactivity of 10 
times bulk earth values, core radius of 550 km, and 
the highest compositional density gradient satisfy 
this constraint as the top of the mixed layer intersects 
the solidus after about 500 Myr. But models with 
smaller bulk radioactivities do not. The models with 
bulk radioactivities of 5 times bulk earth values 
generate melts only at more than one Gyr after core 
formation. Models with small compositional density 
gradients have similar difficulties. For the same bulk 
radioactivity and mantle viscosity, the models with 
0.3%, 0.6% and 0.9% density change per 100 km 
generate melts at 1000 Myr, 636 Myr and 507 Myr, 
respectively. It is clear that the low gradient case 
cannot produce melts in the time required and the 
intermediate gradient case is only marginally accept- 
able. Only the highest gradient case clearly satisfies 
the known age constraints on mare volcanism. The 
same is true for the model with a 475 km core 


radius. This would, however, require higher heating 
rates than the model with a 550 km radius core. In 
models that produce melting at the top of the mixed 
layer in 500-600 Myr, the pressure of melting ts in 
the range of 3-3.5 GPa (see I able 1). 

Any model for the generation of high-li0 2 mare 
basalts from such depths faces a serious obstacle if 
high Ti0 2 mare basalt liquids are neutrally buoyant 
with respect to an olivine-rich mantle in the 2-2.5 
GPa pressure range [4,49]. High Ti0 2 mare basalt 
liquids that are more dense than an olivine and/or 
orthopyroxene residue at higher pressures should 
therefore sink rather than rise toward the surface to 
be erupted. One solution to this dilemma is to sup- 
pose that high Ti0 2 mare liquids do not wet the 
grain boundaries of olivine and/or pyroxene. In fact 
Toramaru and Fujii [50] have presented experimental 
evidence that terrestrial ocean floor basalt liquid 
does not wet the grain boundaries of enstatite in 
mantle peridotite. In dry experiments in the olivine— 
orthopyroxene-plagioclase system, melt was also 
excluded from pyroxene grain triple junctions [51]. 
Because the degree of melt connectivity depends 
critically on the composition of the melt, we cannot 
conclude that this behavior applies to high-Ti0 2 
mare liquids in the lunar mantle. But if mare basalts 
liquids are not interconnected for small degrees of 
melting, melt would initially be trapped along grain 
junctions containing orthopyroxene. Such melt would 
be carried upward by upwelling mantle even though 
the melt is not buoyant relative to the mantle in 
which it resides. 

An alternate solution [3] is that the rate of mantle 
upwelling exceeds the sinking rate of negatively 
buoyant melt. Near solidus melt crytallizes as it 
descends along its adiabat. Only above the neutral 
buoyancy depth can melt segregate into descrete 
bodies in which it becomes effectively isolated from 
the mantle. If this model is correct, the high-Ti0 2 
melts are products of polybaric melting and the 
apparent pressure of multisaturation is the pressure 
below which the melt becomes positively buoyant. 


5, Summary 

This paper has explored the chemical and thermal 
consequences of an internal evolution model ac- 
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counting for the possible role of chemical buoyancy 
and differences in radioactive element content result- 
ing from the crystallization of t he lunar magma 
ocean. Rayleigh-Taylor instability causes the dense 
ilmcnite-rich cumulate layer and underlying Fe-rich 
cumulates that form near the end of magma ocean 
solid ifcation to sink toward the center of the Moon. 
After this overturn, radioacitve heating within the 
ilmenite-rich cumulate core heats the overlying man- 
tle, causing the development of a heated, convec- 
tively mixed layer above the core-mantle boundary 
that eventually becomes thick enough and hot enough 
to melt. This deep high-pressure melting, as required 
by the chemisty of the mare basalts, occurs after a 
reasonable time interval to explain the onset of mare 
basalt volcanism if the content of radioactive ele- 
ments in the core is about 10 times that of the bulk 
Moon (taken to be equivalent to bulk earth values) 
and the chemical density gradients above the core 
are comparable to values that would be expected if a 
minimal amount of mixing occurred during large 
scale overturn of the initially unstably stratified lunar 
interior. 

Regardless of details implied by particular model 
parameters, gravitational overturn driven by the high 
density of magma ocean ilmenite-bearing cumulates 
should concentrate high-Ti0 2 mare basalt sources, 
and probably a significant fraction of radioactive 
heating, toward the center of the Moon. This will 
have important implications for both the thermal 
evolution of the Moon and for mare basalt genesis. 
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